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Abstract 

The background of this work is the problem of reducing the aerodynamic turbulent friction drag, 
which is an important source of energy waste in innumerable technological fields (transportation being 
probably the most important). We develop a theoretical framework aimed at predicting the behaviour 
of existing drag reduction techniques when used at the large values of the Reynolds numbers Re which 
are typical of applications. We focus on one recently proposed and very promising technique, which 
consists in creating at the wall streamwise-travelling waves of spanwise velocity (M. Quadrio, P.Ricco 
& C.Viotti, J. Fluid Mech. 627, 161-178, 2009). 

A perturbation analysis of the Navier-Stokes equations that govern the fluid motion is carried out, 
for the simplest wall-bounded flow geometry, i.e. the plane channel flow. The streamwise base flow 
is perturbed by the spanwise time-varying base flow induced by the travelling waves. An asymptotic 
expansion is then carried out with respect to the velocity amplitude of the travelling wave. The 
analysis, although based on several assumptions, leads to predictions of drag reduction that agree 
well with the measurements available in literature and mostly computed through Direct Numerical 
Simulations (DNS) of the full Navier-Stokes equations. New DNS data are produced on purpose in 
this work to validate our method further. 

The method is then applied to predict the drag-reducing performance of the streamwise-travelling 
waves at increasing Re, where comparison data are not available. The current belief, based on a 
.Re-range of about one decade only above the transitional value, that drag reduction obtained at low 
Re is deemed to decrease as Re is increased is fully confirmed by our results. From a quantitative 
standpoint, however, our outlook based on several decades of increase in Re is much less pessimistic 
than other existing estimates, and motivates further, more accurate studies on the present subject. 

keywords: Drag reduction, channel flow, moving walls, asymptotic expansions 

1 Introduction 

Driven by strong technological interest, the fundamental problem of manipulating turbulent flows is re- 
ceiving more and more attention by the scientific community Significant leaps forward in our physical 
understanding of turbulence and in our ability to simulate turbulent flows, either numerically or experi- 
mentally, have contributed in recent years to raise even further the interest in such topics, together with 
the growing concern with the energetic issue and environmental pollution. 

Perhaps the most difficult problem in this field is the control of turbulent wall flows to the aim of 
reducing the skin-friction drag. Applications (for example air-, water- and ground-based transport, as 
well as duct flows) are wide and the scientific challenge is significant, since wall flows, even in their 
simplest geometry, contain the essence of turbulence, i.e. the near-wall regeneration cycle [51 117) which 
adsorbs energy from the mean motion and redistributes it in an anisotropic way through the mediating 
action of pressure among turbulent fluctuations of the various velocity components. 

A number of strategies to attack the problem exists, ranging from the so-called passive techniques, 
like riblets [3j which modify the geometry of the planar wall in the hope of improving performance, to 
active, feedback-control techniques [5] which build upon linear control theory to devise a control kernel 
capable of driving a large number of distributed microactuators on the basis of an input given by a large 
number of distributed microsensors. The simplicity of passive techniques is counterweighted by their so 
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far limited performance; on the contrary, feedback control is performing well (in numerical simulations), 
but implementation issues are overwhelming. 

The best of both worlds could be found in the intermediate group of active, open-loop techniques. As 
for the passive techniques, they are reasonably simple to implement, not requiring distributed microsen- 
sors and microactuators. As for the feedback-control techniques, their global performance may be good 
enough to produce potential energy savings that are significantly larger that their working cost. One 
recent and promising strategy, especially so owing to its interesting performance in terms of net energy 
savings, has been described for ducts flows and boundary layers, and consists in creating at the wall a 
suitable distribution of spanwise velocity. The spanwise velocity varies sinusoidally in space as well as in 
time, to produce a spanwise-uniform wave that travels along the streamwise direction. This technique, 
introduced by Quadrio and coworkers in 115] and reviewed by Quadrio in |13) , can fully relaminarize the 
flow at low values of the Reynolds number Re, and yields net energy savings of more than 20% for the 
higher Re tested so far. An experimental verification exists pQ, where drag reductions near to 50% were 
measured in a low-i?e turbulent pipe flow. 

As underlined in [T3], one key question that has never been seriously addressed so far, however, is how 
these techniques behave when Re is increased from the low values typical of the available experimental or 
numerical analyses to the high values typical of the applications. Direct Numerical Simulation (DNS) has 
been so far the numerical technique of choice for the minimum amount of flow modelling implied; however, 
it is obvious that, given the tremendous increase of its computing cost with Re, DNS will not allow us 
to get the much needed high-i?e information. Experiments here are stuck to an identical impasse, since 
the not-yet-satisfying development of suitable actuators limits the experiments to similarly low values of 
Re. Two scenarios are possible, and consistent with the limited data available, which evidence a mild 
decrease in performance (identified with the maximum obtainable drag reduction) over the extremely 
modest range of Re explored so far (less than one decade above the subcritical value). One is, of course, 
that performances keeps decreasing as Re increases, thus leading rapidly to a no-benefit situation well 
before the application-level Re are reached. The other suggests that the observed decrease concerns the 
low-i?e regime only: at higher Re performance stabilizes or at least decreases very slowly. This last 
scenario is consistent with the observations put forward in by Iwamoto et al., who studied the effect 
of completely removing the turbulent fluctuations in a thin near- wall layer, whose thickness was kept 
constant in wall units, and observed that the consequent reduction in skin-friction drag, after the initial 
drop at low Re, remains quite high even for high values of Re. 

Early LES-based results are beginning to appear |18) . and seem to present analogous problems, pro- 
viding reasonable results only when the employed spatial resolution becomes comparable to that of DNS. 
It becomes thus obvious that for understanding the high-i?e regime one has nothing left but resorting to 
the numerical solution of the Reynolds-averaged Navier-Stokes (RANS) equations closed with a suitable 
turbulence model. Unfortunately, no model exists to date that is capable of accounting properly for the 
new physics brought about by the wall-based control, inducing such large reductions of turbulent drag. 
The approach recently proposed by Moarref and Jovanovic [TT] resorts to calculations based on RANS 
equations where an eddy viscosity is computed on the basis of spectral information of the non-controlled 
flow obtained by DNS. Though interesting in principle, and applied so far to the oscillating wall only, their 
approach presents the drawback of requiring prior DNS information, and thus precludes its applicability 
to the high-i?e regime of interest here. 

In this paper, we follow an approach similar to [11] and present a predictor model that is based 
on the RANS equations and aims at predicting turbulent drag reduction without requiring high-i?e 
DNS information. Our model is specialized to the case of the streamwise-travelling waves, and how it 
performs in different situations still has to be verified. Moreover, it is a very simple model that makes 
several assumptions, some of which are known to be not entirely justified. It is however the best we 
can presently do by a perturbative approach, and we will demonstrate in the paper that, when properly 
developed by taking advantage of physical insight, the model provides us with useful results that may 
serve as an important guideline for improving our understanding of the problem at hand. The main 
assumption of the model is that the effect of the wall-based travelling waves is confined near the wall. 
This is a very reasonable assumption, as it is shown in [T3] that the waves produce a Generalized Stokes 
Layer (GSL), which is a generalization of the conventional Stokes layer produced by a wall in harmonic 
motion under a still fluid, and that the thickness of such GSL, compared to the distance between the 



2 



channel walls, is extremely small when the waves produce drag reduction. Under this main hypothesis, we 
try and represent the effect of the GSL through a turbulent viscosity that depends on the wall conditions 
through a perturbative hypothesis that is applied both to the velocity field and to the turbulent viscosity 
itself. 

The structure of this paper is as follows, ^describes the perturbative approach, with basic equations 
introduced in 32.11 the zeroth-order problem discussed in 32.21 and the first-order one in 32.31 the main 
quantity related to relative drag reduction is discussed in 32.41 Next, $3] presents results of our predictive 
model and compares them with available DNS information. Then, 2] uses the predictive capabilities of 
the model to extract new information on the behavior of turbulent drag reduction, first in t j4.ll at high 
wavenumbers, and then in 34.21 at higher values of Re. Lastly, t}S]is devoted to a concluding summary. 



2 The perturbative approach 
2.1 Basic equations 

The physical problem under study is characterized by a wall forcing, consisting in a spanwise velocity 
distribution which varies sinusoidally in space and time with a spatial scale L wa u and a time scale 
T wa ii- As a consequence, this distribution travels along the streamwise direction with a speed given by 
L W aii/T wa u- This symmetrical forcing does not introduce any net mean flow in the spanwise direction. 

The definition of the mathematical model under consideration here relies first on the following defi- 
nition for the mean value of a generic function of time f{t): 

i r T 

1 = j I /*, (i) 

where T <C T wa u, i.e. T must be large enough to smooth out the turbulent fluctuations, but it does not 
filter out the wall movement. As we will see in the following, this hypothesis is justified as the travelling 
waves in their most interesting drag-reducing regime are characterized by a relatively long time scale 
compared to the turbulence time scale. 

The basic equation system required here is the standard three-dimensional momentum and continuity 
equations set: 

d t Vi + vjdjVi = --dip + {vdjVi - v[v'^j (2) 
dm = (3) 

The geometry and reference system are shown in figure [I] x, y and z denote the streamwise, wall- 
normal and spanwise coordinate, with u, v and w the corresponding velocity components. We consider 
an indefinite turbulent plane channel flow, for which x and z are homogeneous directions, so that d x and 
d z are null for any velocity statistics. The flow is driven by a constant longitudinal pressure gradient 
d x p — Gp < 0, whilst the wall conditions are u — v = (and w — if no travelling waves are applied). 
The basic equations can be easily made nondimensional on the basis of the channel half-width S and 
the bulk velocity. Here, to ease the comparison with [T3], we choose as a reference velocity scale the 
centerline velocity U of a Poiseuille laminar flow having the same mass flow as the channel under study. 
This velocity scale is related to the bulk velocity by a factor 3/2. 

After the introduction of the standard Boussinesq hypothesis, which expresses the Reynolds stresses 
as proportional to the gradient of the mean velocity field, the coefficient of proportionality being the 
turbulent viscosity: 

2 
3 

the equations ©, © become 



v'^j = -KSij - v T {djVi + diVj) , (4) 
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d t u + ud x u + vdyU = — d x p + -d x n + [y + Ut){9 x u + dyU) + 2(d y VT)d x u + (d y VT)(d y u + d x v) (5) 
P 3 
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d t v + ud x v + vd y v = -d y K + (v + v T ){d x v + dyv) + 2{d v u T )dyV + {d x VT){d v u + d x v) (6) 
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Figure 1: Sketch of the coordinate system and the flow under study. The waves of spanwise velocity (with 
wavelength 27r/fc and oscillating period 2tt/u)) are applied at both walls and travel along the streamwise 
direction with speed c = uj/k that can be cither positive (forward-travelling wave) or negative (backward- 
travelling wave). 



d t w + ud x w 4- vdy-w = (y + i>T)(d x w + d 2 w) + {d x VT)d x w + {d y VT)d y w (7) 

d x u + d y v = . (8) 

The perturbative hypothesis considered here expresses the mean turbulent flow as a basic part (the 
standard channel flow) plus a perturbation induced by the moving walls (a forcing acting on the channel 
flow) and can be written as: 

u = u (y) +eu 1 {x,y,t) +0(e 2 ) (9) 
v = + evx(x,y,t) + 0{e 2 ) (fO) 
w = + swi(x,y,t) + 0(e 2 ) , (fl) 

where the perturbations are due to the boundary condition expressing the presence of the streamwise- 
traveling waves: 

w(x, 0, t) = £ e ^ kx -^ = ee w (12) 

that is u>i(x,0,i) = e* 6 *, scaled by e — w max /U <C 1. Here the physical oscillation is represented by 
the real part of e l6 , and in what follows observable quantities are represented in general by the real 
parts of functions whilst the imaginary parts may contain related and/or redundant information. The 
small parameter e is assumed in this perturbative approach as the typical velocity amplitude of the 
superimposed w-perturbations, made dimensionless with respect to U. Given the typical values of w max 
described in [T5], this assumption is not entirely justified, since w max /U can reach up to unity. As 
customary in asymptotic expansions |19l 17]. however, the small-amplitude assumption is useful to verify 
whether the analysis is able to yield reasonable results and predict important trends. 

The form of condition (fT2"|) suggests to hypothesize the functional form w = e f(y) e l6 + 0(e 2 ) for w, 
and then to search for solutions of the kind 

u = Uo(y) + e g(y) <f> x (6) (13) 
v = + eh(y)4> y (9) (14) 
w = O + ef(y)e 10 (15) 



4 



at the first order in e, where f(y) is normalized to the w-perturbation small amplitude e at the wall, so 
that /(0) = 1. The functions 4>i can be expressed as generalized Fourier series 

= c + Cl e 4e +c 2 c 2ie + ... (16) 



and the simplest model is obtained by truncation at the first order: 

u = u {y) + e g(y) (a + a ie lS ) (17) 
v^0 + eh(y)(b o + b 1 e' e ) (18) 
w = + ef(y)e ie . (19) 

By substitution in the continuity equation (JSJ it can be easily shown that solutions of this kind can exist 
only if bo — 0. Coefficients a\, b\ can be factorized, so the model takes the form 

u = u Q (y)+eg(y)(a + e w ) (20) 
v = + eh(y)be ie (21) 
w = + ef(y)e l9 (22) 

where the functions f,g,h and the constants a,b may depend in general on the parameters k and w. 
Since the forcing modifies the mean velocity gradients, the perturbative hypothesis can be extended in a 
similar way to the turbulent viscosity: 

v T = Mv) + evi(y) ( 23 ) 

where also v\ may depend on the parameters k and lj defining the forcing. 

Equating the terms of order e° and solving the ensuing system leads to the determination of the basic 
flow, i.e. the channel flow without forcing, which will be briefly discussed in £12.21 Solving the problem 
corresponding to higher order in e eventually provides us with the functions /, g and h, so that the drag 
reduction effect can be discussed. This will be the subject of ; 



2.2 Order e° 

Dimensional analysis on equations ©-([7]) may be carried out by considering at first the unperturbed 
channel flow velocity and length scales U, 8 as in figure [T] together with the streamwise length L, such 
that L/J > 1. At the perturbative level, the motion has the new scales e for the velocity amplitude, and 
T wa ii ~ L wa u ~ 1/k for the variations. We make the hypothesis that the wave velocity is of the 
same order of the stream velocity, i.e. L wa u/T wa u ~ U; again, this assumption is justified by considering 
the typical working conditions of the streamwise-traveling waves |15j . Then the growth order analysis 
leads to 

u-equation = 0(U 2 /L) + 0{e)0{U 2 k) (24) 
w-equation = 0(5/L) [0(U 2 / L) + 0{e)0{U 2 k)} (25) 
w-equation = 0(e) 0(U 2 k) (26) 

Thus, at order e° (unperturbed channel flow) the problem is properly described just by the dominant 
0(U 2 /L) terms of eq. (|5|), whilst the higher order equations and terms will enter at the first order of 
approximation. 

After substitution of equations ([2U]) - (j2"5j) in the the e°-terms of eq. © give 

-G + v' o u' o + (v + v o )u% = (27) 

(the symbol ' denotes y-derivatives) where G is the longitudinal pressure gradient term, whilst at the 
same order the w-equation is negligible and the iw-equation is trivially 0=0. Eq. (I2T1) can be integrated 
giving 

A-Gy + {v + v ) u' = (28) 
where the constant A is related to the wall condition, provided that i*r(0) = ^o(0) = 0: 

A = -vu' a (0), (29) 
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i.e. A is proportional to the basic wall stress. For the present purpose, it is sufficient to consider a single 
wall, setting there the origin of the y axis. The e° problem is then closed by the boundary conditions 
u = v = w = Q&t the wall and by the knowledge of one of the functions Vo(y)> u o(y)- The solution 
can be obtained by standard procedures, namely i) an hypothesis on the shape of vo(y) gives uo(y) by 
solving eq. (|28l) with uq(0) = 0; or ii) the knowledge of Uo(y), from numerical simulations, experimental 
fits and/or special functions approximations, that leads to Vo(y) through equations (|2"8" |) .(|2"9" ]l : 

Gy + is[u o '(0)-u Q '] 

— ; . (OU) 

u' 

The resulting channel flow (unperturbed, e = 0) will take the form 

u = u (y) (31) 
v = (32) 
w = (33) 

where uo(y) represents the progressive velocity increase from the wall, loosing validity near the channel 
center where the derivative of the real velocity profile vanishes, a behaviour that cannot be reproduced 
here because eq. (f2"5|) degenerates for u' (y) — 0. This is not a limit in this work, since the center region of 
the channel is not playing a key role in such wall-based drag reduction techniques that aim at modifying 
the near- wall turbulence regeneration cycle [5] , and the drag reduction effects are supposed to depend on 
the interaction between the moving wall and these flow structures |13j . 

2.3 Order e 

The dimensional analysis tells that the dominant momentum eqs at the e order are the u— and in- 
equations, that can be obtained looking for e-terms in the eqs ([S]),...© after substitution of eqs (|2T)|) . . . . (|2U|) . 
The equation system to consider now consists of the w-equation 

(l + ae~ ie ) (v + v ) g" + (l + ae~ ie ) v' g' + (i u - i k u - k 2 v - k 2 v ) g + 

+ (ikv' -<u' ) bh + e-^iv'^ + vxu'^) =0, (34) 

the ui-equation 

{v + vq) f' + vof- (k 2 v + k 2 is Q + iku -iuj) f = (35) 
and the continuity equation 

9 = i^' (36) 

(the effects of k in the denominator will be discussed below). 

Equation (|34p involves the unknown functions g, h and v\. Substituting (|36[) in (|34[) the following 
equation for h is obtained: 

c~ 10 [iab [v + vq) h'" + i ab v' h" + fc i>[ u' + fc v\ u'q ] + 

+i (y + v ) ti" + i v' Q h" + (k u - u - i k 2 v - i k 2 v ) ti + (i k 2 v' Q - k u' ) h = . (37) 

This equation is of the kind e~ 10 [C\(y)h(y) + M.\(y)) + [£,2{y)h{y)} = 0, where C, M. are linear operators, 
and proper solutions exist only if both parts vanish. The first part C\h + Ai± =0 is 

i a b (v + uo) h!" + ia b v' Q h" + k v[ u' + k V\ u'q = 0, (38) 
and can be integrated to give 

B + k^u'a + 'iabiv + v^h" = (39) 
(B is a constant). This relation is satisfied by 

Vx = _£±i^+i^! ! (40) 
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and the value of B can be obtained through the condition fi(0) = 0. 
The second part L^h = of eq. (j37|) . namely 



i (v + uq) b!" + i v'q h" + (k i*o — ui — i k 2 v — i fc 2 ^o) /V + (i & I'o — ^ u o) ^ — 0, (41) 

can be solved for h with boundary conditions h(0) = h'(Q) — 0, i.e. iti = v\ = at the wall, together 
with a third condition. This gives 

h = h{y;k,LU,s) (42) 

i.e. the function h depends on y, on the parameters k and u>, and on a integration constant s. 

The w-equation (f3"5)l may be solved for / with condition /(0) = 1, i.e. w = ewi — ee 10 at the wall 
after eq. (TTTj) and (TT5)) . together with a second condition. 

We note here that, for a laminar flow, the equations for u and w become simpler. In particular the 
laminar u-equation (|37|) reads 



e 10 i a b v h'" + i v h'" + [(k uq — to) — i k 2 v\ h! — ku' h = (43) 

i.e. takes the form e~ 10 [£i (y)h(y)] + [£ 2 (y)M^)] = 0- Hence, the first part C\{y)h{y) = must be 
satisfied by a trivial condition, i.e. a — and/or 6 = 0; this condition will be identified in the next 
section. The second part £2{y)h(y) = could be solved for h as in the turbulent case, with the same 
boundary conditions, but in the next section it will be shown that the solution of this equation is not 
necessary in the laminar case. 

For a laminar flow, even the w-equation (|35p becomes simpler: 



vf- [k 2 u + x{ku Q -io)\ f = (44) 

and admits decaying solutions in terms of confluent hypergeometric functions. For a very short y-domain, 
where the parabolic expression of uo (the Poiseuille velocity profile) can be linearized, the solutions take 
the form of Airy functions, in agreement with the analytical solution for the laminar GSL determined in 

El- 

2.4 Drag reduction 

The total wall stress is, after equation (1201) . 

T = p\v + v T {Q)\ (d y u) y= o =pis[u' o {0)+eg'(0)(a + e w )] , (45) 
and in nondimensional (macroscopic) form: 

The drag-reducing effect of the travelling waves can be estimated by calculating the mean value c/ >m of 
the friction coefficient c/ over a large time T g such that T g 3> T wa u, this quantity is 

c/,m = p[«o(0)+«ff'(0)]- ( 4 7) 
The unperturbed wall stress is 

c fo = ^u' (0), (48) 
so that the relative drag reduction, after equation (|36l) . is: 
r.R i c f> m £a 9'(0) ieabti'jO) 

UK=1 = T7n\~ = 1 T7n\ ' ^> 

c f0 u' (0) ku' o (0) 

which means 

DR(fc, u>) oc ^ h"(0; k, uj, s) (50) 
k 
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where C — ab. The relation (|50p allows a map of the relative drag reduction to be drawn as a function 
of k,ui and of an integration constant s, which has the meaning of an accessory condition (obviously, 
the dependence of drag reduction upon the amplitude of the spanwise forcing cannot be determined by 
this perturbative approach). The factor C cannot be determined at this level of approximation, but the 
required regularity as k — > of g'(y) and DR(fc, u>), both proportional to b/k after eq. ([55]) and ([501) . 
suggests to assume 

b(k) ~ k n with n > 1 as k -> 0, (51) 

an hypothesis which can be formally introduced in eq. (|20p . a priori. 

It is important to note here that the physical meaning of DR is expected to be found in its real part, 
which depends both on the real and imaginary parts of h since the factor C is in general a complex 
quantity. The real part of eq. (|50p . writing C = C r + iCi and h = h r becomes 

DR r (/c, w) oc - [C r h"(0; k, u>, s) - Ci h'!(Q\ k, w, s) 1 (52) 
A; 

where we have assumed C = a 6 ~ fc™ with n > 1 as — >• as stated by condition (fBTj) . 

For a laminar flow, a and/or 6 must vanish as explained in 32.31 so that equation ((49]) gives DR = 0. 
Again, this result is consistent with the results described in [14], where the laminar Poiseuille flow is 
demonstrated not to be affected by the GSL originated by the travelling waves, since the streamwise 
parabolic velocity profile and the transverse GSL profile are entirely decoupled. In our model, the basic 
equations ©-((HJ) for a laminar flow are no longer coupled through the turbulent viscosity vt, and this 
suggests reconsidering the basic perturbative assumptions (|20[) -([22 p . checking their appearance when 
a = 0or6 = 0ora = 6 = 0. It turns out that the right closure of the problem for the laminar case is 
6 = 0, because in this way the assumptions (I2TJ1) . . . (|2"2"|) take the form 

u = u o(y) (53) 
v = (54) 
w = ef(y)e te (55) 

(6 = in the continuity equation (1361) gives g(y) = 0) and represent a flow where the x and y-momentum 
equations are independent from the w velocity component, whilst this component can be calculated by 
solving equation (j44]) . that couples / with uq. 



3 Results 

In order to obtain h"(0) and thus the map of the drag reduction DR as a function of the parameters k 
and u, eq. (|4ip must be completely characterized by specifying the functions uq and isq, then solved with 
the relevant boundary conditions on a suitable domain < y < y m , i.e. a range of distances from the 
wall to a position where the drag reduction effects of the moving wall become negligible. The details of 
this procedure are presented in the following sections. 

The solution procedure is expected to lead to a DR(fc, ui) map, that requires first to be compared with 
known results for validation. Quadrio et al. in [15] provide an ample dataset created with a comprehensive 
numerical study. They carried out state-of-the-art Direct Numerical Simulations (in terms of computer 
code and numerical algorithms, and - most importantly - in terms of discretization parameters and 
solution procedure) to determine the dependence of drag reduction DR upon k and u) at Re — 4760. 
Fig. [5] reports their results, and can be considered as a reference map for comparison. This map clearly 
exhibits an oblique region near the vertical axis where DR reaches its maximum value (referred to as 
hill in the following) and a triangular region where DR has a negative minimum, being a region of drag 
increase (valley). These regions are separated by an oblique line where DR = 0, called the neutral line. 

3.1 Problem closures and general properties 

Once Mo and isq are known, the problem at hand is given by eq. (|4ip together with its boundary conditions. 
The most general way to impose these conditions is based on the consideration that the physical meaning 
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Figure 2: Map of drag reduction (percentage) at low Re in the k-uj plane, adapted from 15 . The 
thick contour line marks the neutral level (no drag variation). Friction drag increases within the darker 
triangular region. 



is contained in the real part of h and g, so the problem may be written in the form 

C 2 h = (56) 
3?{/i(0)} = (57) 

Hh'm = -£»{$(0)} = (58) 
S*{fc'(Wm)} = or 5?{% m )} = (59) 

where conditions (fST)) and (fS"8")) state that the real parts of h and g vanish at the wall, i.e. u\ = v\ = 0, 
whilst condition (|59[) states that g ox h {u\ or v\) vanishes at the far boundary y m , representing 2 possible 
closures of the system. A first study of this system can be made substituting power series for uq and vq 
of the kind 

JV N 

u a = Y J *ny n + o(y N+i ) ] ^o = E 6 »y" + °^ Ar+1 ) ( 6 °) 

n=Q n=Q 

and according to these approximations, a similar series for the unknown function h: 

N 

h = Y / c n y n + 0(y N+1 ). (61) 

n=0 

The resulting set of equations gives non trivial solutions (where any c„ is non zero) only if h and g at 
the wall vanish completely (real and imaginary part), so that the problem under study becomes 

C 2 h = (62) 

h(0) = (63) 

h'(0) = (64) 

%{ti(y m )} = or 3?{% m )}=0. (65) 
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The two closures expressed by condition (|65J) become in explicit form 

h'(y m ) =s or h(y m ) = is (66) 

where s is a real number with the meaning of an integration constant. 

Several properties of the function h can be revealed by the power series method (see appendix lA"j) . 
The most important are outlined in what follows for the problem closed by condition h'(y m ) = s; it is 
easy to see that these properties remain unchanged even with the second closure, with the exception that 
the roles of real and imaginary parts h r and hi are exchanged, as well as their signs. The power series 
method gives general results of the kind 

h"(0;k,uj, S )^ — — (67) 

P(k,uj;y m ) 

where P{k,ui;y m ), having y m as a parameter, is a k, w-polynomial whose order depends on the chosen 
order for the function h. The parameter s turns out as a factor owing to the homogeneity of eq. (|o12)) . and 
only affects the amplitude of h"(0; k, uj, s), that remains unknown anyway because of the perturbative 
approach. From ((67)) it is easy to see that h"(0; k, uj, s) is regular and tends to as k — > ±oo and ui — > ±oo. 
This property has just a mathematical meaning, since under this limit the wavelength and period of the 
wall forcing L wa u and T wa u tend to 0, which is a violation of the initial hypothesis that the forcing takes 
place on scales that are large with respect to the turbulent spatial and temporal scales. 

After defining the domain < y < y m , it can be seen that the use of a truncated power series for 
h of order > 5 is sufficient to reproduce the basic features of the original, DNS-computed map of DR, 
leading to imaginary parts h" that exhibit a oblique hill and a parallel valley separated by a neutral 
line. Actually, once h is known, a reconstruction of the DR map may be attempted using eq. (I52[) . that 
becomes here 

DR r (fc, uj) oc - [C r ti'(0; k, uj) - CU"(0; k, uj)} 
k 

where the constant s is no longer highlighted thanks to its factorization. Another general property tells 
that in the origin of the k — uj plane h'( — whilst /i" ^ 0. This suggests to assume C r = since at 
(k,uj) = (0,0) the wall motion does not exist and in this case the drag reduction must reduce to zero. 
The dependence of C on k and uj can be expanded in a standard power series of the kind 

C = Co + Ci k + C 2 uj + C* 3 k 2 + C 4 k uj + C 5 uj 2 + ■ ■ ■ 

that becomes, owing to hypothesis (|51l) about regularity in the origin, 

C = k{C 1 +C 3 k + C 4 uj + ■■■). (68) 

The practical use of this expansion requires a truncation, and it can be seen that a linearized form of 
C(fc, uj), together with a truncated series of order > 5 for h, are sufficient to give an approximated version 
of the DR map comparable to the reference one. A similar method will be applied in the following sections 
starting from functions h numerically calculated, and leading to results of better accuracy. 

3.2 Numerical solution 

The search for a numerical expression of the function h(y) starts by considering the original equation 
(|¥Tj) together with the accessory conditions, i.e. the problem (IB^I) . . . ([55)1 . The functions uq and Vq 
which characterize equation (|4ip are defined by the method (ii) outlined in H2.2[ starting with a formula 
for uq expressed in terms of special functions that accurately approximate the mean velocity profile, 
which is linear very near to the wall (viscous layer), then gradually changes through the buffer layer to 
a logarithmic behaviour which is typical of the so-called logarthmic layer. 

With reference to the accessory conditions, it must be noted that the numerical approach allows an 
easier closure for the original problem, which can be imposed at the wall, in the following way: 

C 2 h = (69) 

h(0) = (70) 

h'(0) = (71) 

h"(0) = S, (72) 
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where S is the integration constant necessary to the closure. Here S is a constant with respect to y, but 
it may depend on the other parameters, i.e. S = S(k,u). 

Then, it is useful to take advantage of homogeneity of eq. C 2 h = 0, which allows to factor out the 
integration constant S. A new normalized function h can be introduced: 

h (y; k, u, S) = S h(y; k, lj). (73) 

This leads to the new system: 

C 2 h = (74) 

h(0) = (75) 

fc'(0) = (76) 

h"(0) = 1. (77) 

Finally, after determining the numerical solution h, it is possible to set S in such a way as to make the 
system (|74 p -(|77 p equivalent to the original one (|62 |) -(|65 p . To do that, S can be defined in two ways: the 
first one is 

S =T ^— so that h'(y m )=Sh'(y m ) = l (78) 
h'{y m ) 

i.e. ^s{h'(y m )} = 0, which is the first closure of the original system already used for power series solutions. 
The second possible definition for S is 

S = y^— so that h(y m ) = Sh(y m ) = i (79) 
h(y m ) 

i.e. $t{h(y m )} = 0, which is the second closure of the original system. It is easy to show that definitions 
of S of the kind S = m/h'(y m ) an d S = \ m/h(y m ) where m is a number are equivalent to the proposed 
ones, since a constant m can always be factorized. 



3.3 Drag reduction 

The numerical solution of the equation for h(y) can be obtained by standard methods (Gear/ Adams). 
The closure condition of the first or second type, respectively, where ^{h'(y m )} or 3?{/i(y m )} vanish at 
a given distance y m from the wall, can now be formally imposed at any y m in the y domain from to 
positions far from the wall. 

As a first step, the problem at the reference Reynolds number Re — 4760 is considered with the first 
closure condition. The function under study is h"(0;k,uj,S) = h"(0; fc, lj, y m ) because the integration 
constant S depends on y m directly, after formula (1781) . An outline of the results for various values of y rn is 
shown in fig. [3] In general, as already revealed by the power series study mentioned in M3 - H the numerical 
solutions too have a real part /i" = 3?{/i"(0; k, u, y m )} that does not tend to zero as (k, u) — > (0, 0) whilst 
the imaginary part does. Again, drag reduction appears to be better represented by the imaginary part 
h'( = 3{/i"(0; k, uj, y m )}, a surface always characterized by a main hill-valley system having minor lateral 
undulations side by side (it can be verified that it tends to as k — > ±oo and lu — > ±oo). This can be 
formally expressed setting C r — in the equation (|5U| . 

At this point, determining the correct value for the parameter y m still is an open problem, that can 
be solved by searching for the y m value such that the neutral line between the hill and the valley in the 
k, lu plane is as close as possible to the one in the DNS reference DR map. This happens for quite a small 
value of y^, namely = 27, confirming that the effects of the forcing are important only in a thin zone 
adjacent to the moving wall, approximately the size of the buffer region of the mean velocity profile. 

The choice of the particular value for y m value is not critical: we have verified that the following 
results are almost insensitive to relatively small changes of y m , as long as y m does not approach the 
central region of the channel, where they start changing dramatically. Furthermore, it can be seen that a 
study based on the second closure condition 3?{/i(y m )} = leads to analogous conclusions. The h" map 
with the first closure for ?/+ = 27 at the reference Reynolds number is shown in fig. 2) 
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-Im{h"(0)} 



-Re{h"(0)} 
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Figure 3: Maps of the imaginary (left) and real (right) parts of the function h"(0; k, oj, y m ), computed 
for various values of the parameter y m , and with the closure condition of the first kind at Re — 4760. 
The neutral line, as all the points where h" = 0, is evidenced in the maps of h" (0; k, co, y m ) as a thick 
solid line. 
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Figure 4: Map of h'((y\k,uj), at the optimal y m 
and at the reference Reynolds number Re — 4760. 



Figure 5: DR r map approximated by formula 
— (C^ik + C^uj) /i"(0; fc, uj) where the coefficients 
C^i and are optimized on the ku domain un- 
der study at the reference Reynolds number. 



The reconstruction of the DR map can now be attempted by using equation (1521) . i.e. 

DR r (fc, uj) oc - [C r h"(0;k, uj) - Q/i-'fOjfc, uj)} 
k 

(the dependence on y m , which is now a constant, is no longer highlighted). As explained above, the 
simplest representation of DR can be obtained by setting C r = 0, whilst Ci(k,uj) must be determined in 
such a way that 

DR r (fc, uj) = —-Cih'-(0;k, uj) 
k 

is as close as possible to the reference map. This can be done by expanding Ci in powers of k and uj after 
equation (|68[) . truncating the series at the minimum possible order and imposing the necessary properties 
in such a way as to determine the unknown constants, i.e. imposing regularity, DR = at (k, uj) = (0, 0) 
and symmetry (k,uj) -H- (— fc, — uj) since the physics of the traveling wave remains the same in both cases. 
The simplest truncation which retains the fc, w-dependency is obtained here by setting C n = for n > 5, 
what gives DR r oc —Cih"(0; k, uj)/k — — (Cu + C^fe + C^uj) h"(0; fc, uj), that vanishes in the origin. The 
search for the best coefficient (Cu + C^k + C^uj) fitting the reference DNS map shows that Cu must 
vanish because of the inherent symmetry of the DR function for (fc, uj) (— k, —uj) . In this way it is 
obtained 

DR r (fc, uj) = -{C 3l k + C 4l uj)h^(0;k, uj). 

By using the function h" calculated as in fig. |4j the optimal values of the coefficients on the k, uj domain 
under test turn out to be C3.; — 0.22, — —0.24. The resulting map is shown in figure [5j and is in good 
agreement with the DNS map, the main characteristics of which are qualitatively reproduced: namely, 
the two-lobed side-by-side pattern with drag reduction and drag increase, and the double local maximum 
of drag reduction on the k = axis, corresponding to the temporal oscillating wall emerge clearly from 
the figure. From a quantitative viewpoint, the neutral line has a slope which is not far from the true one, 
and the lobes of drag reduction and drag increase have position and extent which are similar to what can 
be observed in the DNS map. 

The C(k,uj) function involved in figure[5]is a compact, linearized expression, but nothing prevents us 
from considering better representations for this factor. Relaxing the assumption C' r — 0, while keeping of 
course the condition DR(0, 0) ~ and searching for higher order approximations, a representation of the 
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Figure 6: DR r map approximated by formula 
i [ C r <(0; k, u) - Ci h'f(0; k, u) } where the coef- 
ficients Ci and C r are expressed as 2nd order k, ui- 
polynomials ratios after formula (|80D and opti- 
mized on the k — uj domain under study at the 
reference Reynolds number. 



Figure 7: DR r values from DNS (dots) and present 
model (continuous lines) based on the same coeffi- 
cient C used in figure [5] at k = 1 (top) and k = 2 
(bottom). 



DR surface very close to the reference one can be obtained, as reported in fig. [5] In this case, C(k,u>) 
still can be expressed as a truncated power series of the kind as before, but a better representation 
can be probably written in terms of rational functions: 

c = k Cox + C a2 k + C a3 uj + C a4 k 2 + C a5 kuj + C a6 uj 2 + ■ ■ ■ 
Cbi +C b2 k + C b3 uj + C M k 2 + C b5 kuj + C b6 uj 2 H ' 

and C r is expressed in a similar way. The functions Ci , C r used in the DR map of fig. [5] are expressed 
according to formula ([50]) in terms of ratios of simple 2nd-order polynomials in k and ui. 

A more quantitative assessment of the results obtained so far can be had by looking at sections of the 
map taken at fixed values of k: for example, fig. [7] shows two sections of fig. |H1 where it can be seen that 
there is good agreement in the positions of maxima and minima between DNS results and the prediction 
of the present model, the only remarkable difference being in the lateral decay rate of the DR values. 



4 Using the model 

Our model has been developed, and it has been verified that it is able to produce data which are broadly in 
line with the available DNS information. As in a typical asymptotic formulation, a part of the information 
required to solve the problem has been introduced as external data from the underlying physics, namely 
the domain width y m and the related third integration constant s or S. It is worth to recall here that 
y m is the range of distances from the wall to a position where the drag reduction effects of the moving 
wall become negligible, whereas s or S have the meaning of an accessory condition in addition to the 
standard no-slip condition at the wall. Another effect of the asymptotic approach is the appearance 
of the multiplicative function C(k,uj)/k in the equation ([50)1 for the drag reduction, a function that 
cannot be intrinsically determined at the first level of approximation on the small nondimensional forcing 
amplitude e. Even if the basic properties of the DR map can be represented by the h"{k,uS) surface, the 
function C(k,uj) was shown to be useful for a "calibration" of the model on the available information, 
approximating it as a first choice in linearized form. 

It is thus of interest to use the model to predict the behavior of drag reduction beyond the available 
DNS information, which is limited to low values of both the wavenumber k of the traveling waves, and of 
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Figure 8: DR r values from DNS (dots) and present model (continuous lines) based on the same coefficient 
C used in figure [5] at k = 10, 20 and 30 from left to right. 

the Reynolds number Re of the longitudinal flow. In the following we will use our model to investigate 
the trend of variation of turbulent drag reduction when k is increased, and - more interestingly - when 
Re is increased. This is the ultimate aim of the present paper. 

4.1 Drag reduction at higher values of k and uj 

A first attempt at using the predictive capabilities of the present model is devoted to investigating 
the behaviour of turbulent drag reduction at high values of k. In the literature, at the reference value 
Re = 4760 there are no such data available for k > 5 and |w| > 3. This first predictive step is intermediate 
and serves us the purpose of testing the model on a region of the parameter space for which it has not been 
calibrated, but for which DNS data can still be produced to evaluate it. We have then purposely run 18 
additional DNS of turbulent plane channel flow at Re = 4760, similar to what has been done in [15] . The 
DNS code, described in |10) , uses Fourier discretization in the homogeneous wall-parallel directions and 
compact, IV-order explicit finite difference schemes in the wall- normal direction, and employs a partially 
implicit time integration method. The computational domain is identical to that employed in the original 
study: L x — 6irh, L y = 2h and L z = 3irh. The spatial resolution is N x = 320 and N z — 320 Fourier 
modes, and N y = 161 points. The main addition to the available data is the range of wavenumbers and 
frequencies that are investigated, which spans the interval 10 < k < 30 and 3 < < 15. 

Owing to the relatively limited number of new data points available, it is difficult to build a reliable 
two-dimensional map. Figure [5] thus compares DNS data (dots) and predictions from the present model 
(lines) at constant k as a function of lo. Again, the positions of maxima and minima are in good agreement 
between DNS data and present model. Overall, both results indicate that the hill- valley structure remains 
elongated and has a rather slow decay in the k — uj plane. 

4.2 Drag reduction at higher values of the Reynolds number 

Applying the present model to investigate the turbulent drag reduction for high values of the Reynolds 
numbers is obviously the main thrust of the present study. Unfortunately, this is prevented by the 
problem of determining the coefficient C at these Reynolds numbers. However, another interesting test 
can be carried out, i.e. computing the function h"(0;k,oj) for larger and larger Re. This can be done 
rather easily by for example keeping the first closure condition at the optimal value y+ = 27. In this 
section only the maps for the imaginary part of h" will be considered, but it can be shown that /i" has a 
similar scaling. 

Figure [9] shows the h'( maps computed with our model at 3 values of the Reynolds number, i.e. the 
reference one Reo — 4760 (top), and the values 10i?eo (middle) and 100i?eo (bottom). Focusing first on 
the left column, the top plot is an extension of fig. H]to the wider domain — 30 < u> < 30, < k < 50 
and has been already discussed. As Re increases, the general shape of the map over the same domain 
seems to change only weakly, if exception is made for an evident increase of its size. This means that the 
characteristic time and space scales shrink as the Reynolds number increases. 

Since we are concerned with near-wall turbulence, it is reasonable to suppose for the phenomenon 
to obey an inner, viscous scaling. This is verified in the right column of figure [9J where the same maps 
are translated in wall units and are shown to keep their general size across a 100-times increase in 
Reynolds number. In particular, the right column is obtained by choosing a domain —0.38 < oj + < 0.38, 
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Figure 9: Maps of h"(0; k,ui) with y+ = 27 at Reynolds numbers Re = 4760, lORe and 100i?e . Left: 
outer units, plots over the domain —30 < lj < 30, < k < 50. Right: inner units, plots over the domain 
-0.38 < uj+ < 0.38, < k+ < 0.026 (see the text). 
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Figure 10: Normalized 
amplitude factors vs. 
Reynolds number. 
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Figure 11: Computed decrease of maximum drag reduction as reported in 
table [T0l (dots). The lines are two fits for these data: dashed line is after 
eq. (|8ip. and solid line is after eq. (|8"2"|). 



< k + < 0.026 in wall units, obtained from the conversion at Re = Reg of the original one used in figures 
[5] to [HI that was —3 < ui < 3, < fc < 5 in outer units. Then, all the maps of the left column are plotted 
again over this new domain in wall units, for all the values of Re, evidencing in this way the morphologic 
analogies. Clearly, the inner units maps are not identical, so that a strict Reynolds-invariance property 
can be excluded; but a property of slow Reynolds-scaling could be revealed by the present analysis. 

Of course this perturbative analysis cannot give any indication on the absolute levels of drag reduction 
that can be achieved. Of particular interest, however, is to investigate how the best performance at low 
Re degrades at higher Re, and how the values of k and tu defining the most performing traveling wave 
are observed to change while Re increases. In order to compare results obtained at different Reynolds 
numbers, a proper amplitude factor must be defined for the h'( maps: since DR(fc, uj) oc C h"(0; k, u)/k 
after equation ([30]). here F(Re) = h" max /k^ max is selected as a factor for the present purpose (ki, max is 
the maximum locus), i.e. a quantity that should give information about the trend of DR as function of Re. 
The results are shown in table [TU] and figure [TTJ where the values of F(Re) normalized on the reference 
value obtained at Reo = 4760 are presented: actually, the trend exhibits a slow Reynolds-scaling, and 
raises a new question, about the existence of an asymptotic limit for large Reynolds numbers. 

Figure [UJ shows also two possible fits for the trend of the calculated points. The first fit, plotted with 
a dashed line in the figure, as function of Re T is: 



where a = 1.31, /3 = —0.0433. The second fit, plotted with a continuous line in the figure, as function of 
Re T is: 



where 7 = 0.682, a = 0.945, f3 = —0.171. Judging from the present data, the second fit appears to be 
better than the first, but it requires the determination of 3 parameters instead of 2, and implies the 
existence of a finite, non- vanishing asymptotic limit. 

The present result confirms, first of all, the expected decrease of the performance of the travelling waves 
with increasing Reynolds number. This is of course a very reasonable effect from the physical standpoint, 
since it reflects the gradual decrease of importance of the viscous near-wall cycle in the overall turbulent 
flow, which is more and more influenced by the large-scale events originating in the outer (logarithmic) 
region [5]. Although no information is available for the travelling waves at Re T > 400, we can use the 
information available for the spanwise-oscillating wall, which is a simpler technique driven by a similar 
physics and has been investigated up to Re T = 1000. Touber & Leschziner [JJ5] for example report LES 



G{Re T ) = F{Re)/F{Re ) = a Re? 



(81) 



G{Re T ) = F{Re) /F(Re ) =-f + aRe^ 



(82) 
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and DNS information to substantiate the claim that maximum drag reduction decreases with a power law 
~ .Re" 2 . A similar result was proposed earlier by [2], although on the basis of a dataset computed for 
100 < Re T < 400. These results are consistent with our findings, although our rate of decay, computed 
on the entire dataset available, seems to be slower, with an exponent of -0.043 only. This rate, computed 
on a smaller dataset for Re < IORbq, turns out to be only slightly faster, with an exponent of -0.051. 
We notice, moreover, that our data lend clearly better support to a decay law of the type (JS2J) - This 
kind of fit would be also more in line with the claim by Iwamoto [4! that the decay of the maximum drag 
reduction is a low-Reynolds effect, or at least that at high Re the performance decay becomes very slow. 

By its very nature the present approach, being based on low-i?e information, is unable to account 
for a possible major change with Re of the layout of the drag reduction map in the parameter space. 
That said, the locus (u>, k)i imax of the maximum values of h'[ max considered above seems to be mildly 
dependent upon the Reynolds number. Actually, in the range Reo < Re < 200i?eo the pulsation varies 
in the range < u)f max < 4 • 10~ 5 , which is in the order of the accuracy of the ODE solver employed, 
whereas the wavenumber shows a very small decrease, namely 0.0089 > k+ max > 0.0081. 

One last comment is in order with respect to the net energy savings brought about by the travelling 
waves. Drag reduction contributes to energy savings through reduction of the pumping power, but 
the energetic cost of control must be considered too. There is information available [TB] only for the 
oscillating-wall case, where the power required to oscillate the wall in a plane channel flow is reported 
to slightly decrease with Re T , being proportional to i?e~ 0136 . This is simply the effect of assuming 
Re = 15.43i?e^. 136 as in [12], since in outer units the required power is constant. It follows that the small 
decrease of drag reduction implies a corresponding small decrease in the net savings. This conclusion is 
expected to hold for the travelling waves too, on the basis of the validity of the laminar GSL solution to 
describe the forcing-induced laminar oscillating transverse velocity profile. 

5 Conclusions 

The present paper has introduced an innovative perturbation method aimed at predicting the turbulent 
drag reduction characteristics of the streamwise-travelling waves of spanwise wall velocity for the geometry 
of a plane channel flow. 

The method, based on a perturbative approach, considers the Reynolds-averaged Navier-Stokes equa- 
tions in the simple geometry of an indefinite plane channel flow. In addition to the usual streamwise 
mean flow, a spanwise oscillating flow with zero mean value induced by the wall forcing is present. A 
perturbation expansion is applied to the mean streamwise and spanwise flow, as well as to the turbulent 
viscosity profile. Under the fundamental assumption that the spanwise flow is small with respect to the 
streamwise flow, the perturbative problem expressed at zero-th order corresponds to the standard tur- 
bulent channel flow, whereas at first-order one obtains a new equation set, which contains the spanwise 
flow and its interaction with the streamwise flow. Our interest is focused on the wall derivative of the 
streamwise velocity profile, a quantity that is easily related to the turbulent friction. The modification 
of this derivative as an effect of the wall forcing is then computed in the asymptotic frame. 

The procedure is first used to reproduce the available information on the turbulent drag reduction 
brought about by the traveling waves. This information describes how drag reduction depends on the 
wavenumber and temporal frequency of the waves, and was previously computed by Direct Numerical 
Simulation of the full Navier-Stokes equations. In the present model, a function appearing as a factor 
in the DR formula was expressed in two possible ways, at first as a linearized expression, then as a more 
accurate formula, in both cases this function was used to improve the quality of results, and the relevant 
coefficients have been obtained by fitting the reference DNS map on a reference k, u> domain. 

Once the procedure has been wholly established, including its parameters, it is used to produce new 
information. It is first tested on a new dataset, purposely produced by DNS during the present work, 
designed to explore the drag reduction induced by the waves at low values of Re but high values of k 
and u>. By keeping the fit coefficients computed on the smaller k,u> domain of reference, our method is 
shown to yield a good prediction. 

Lastly, the main question that motivates the present method is addressed, and the variation of the 
maximum drag reduction brought about by the travelling waves when the value of Re increases is pre- 
dicted. This prediction instead has been obtained avoiding the use of fitted coefficients. Compared to the 
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baseline value Reo, the available literature information is limited to 2i?eo for the streamwise-travclling 
waves, and up to 5i?eo for the spanwise-oscillating wall. Here we can easily reach more than 100i?eo, 
so that the two-decades span allows to draw some clear conclusions about a true high-i?e trend. Our 
findings are, first of all, that maximum drag reduction decreases, as expected on the basis of physical con- 
siderations, since the near-wall layer becomes less and less influential on the whole turbulence dynamics 
as Re grows. The important result, however, is that our data suggest a significantly milder decay when 
compared to the available predictions of the kind ~ Re~ const - , which are however limited by the very 
small extent by which Re has been increased. Instead of the commonly reported decay ~ Re~ , we 
suggest a much slower decay ~ Re~ a o . When an asymptotic value for the drag reduction at large Re is 
allowed by the chosen fit, our data show a better agreement, thus implying that drag reduction does not 
decrease below a certain threshold as the value of Re increases. 

The predictive procedure described in this paper undoubtedly contains a number of significant and 
critical assumptions. Just to name the most critical one, we remind the reader of the Boussinesq hy- 
pothesis; of the assumption that the amplitude of the transversal velocity waves is small compared to 
the longitudinal velocity scale; of the assumption that (space and time) scale separation exists between 
the waves and the turbulent flow. That said, however, the scenario that emerges from using the pro- 
cedure is generally consistent, as far as predictions can be critically evaluated against available data. 
When the procedure is used to make true predictions in a regime where no information is available, the 
conveyed message is reasonable and points to a high-i?e scenario that is less pessimistic compared to 
what is generally thought of on the basis of the available information. Without counting too much on 
the exact validity of the predicted scenario, we would like to offer this more optimistic high-i?e view 
to the flow control community, as a further motivation to intensify the efforts towards understanding 
what really happens to wall-based turbulent drag reduction techniques when the values of the Reynolds 
number become significantly high and reach application-level. 
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A Drag reduction from power series method 

The simplest attempt to solve problem (|6"2"j) . . . (|6"5j) without loss of physical meaning makes use of series 
of the kind (|60|) for uo, vq, truncated at the 3rd order in y. The relevant polynomial coefficients a n and 
b n can be fit on the numerous literature data. This level of approximation permits to write h"(0) in a 
closed form of acceptable size, but it is limited to low values of y in the order of < y + * 50, since it 
may become inaccurate over a large domain. As regards the unknonwn function h, it is easy to see that 
the solutions of the problems based on the power series (|6"Tj) for h truncated at the iV-th order leads to 
trivial results until N < 4 with either closure condition, ^s{h'(y m )} — or 3t{h(y m )} = 0. 

The first interesting case is obtained at order 5 in y: here the problem closed by condition *<s{h'(y m )} — 
0, formally expressed as h'(y m ) = s where s is a real number, gives 

24 s is 3 1 

h"(0; k, u>, s) = — (83) 

y m a + 13 w + 7i k + 72 k A 

with 

a = 6 b\ yl - (8 b\ y 2 m + 12 b x b 2 y 3 m ) v + (12 b x y m + 8b 2 y 2 m + 6 b 3 yf n ) u 2 - 24 v 3 (84) 
/3 = My 2 m v 2 - MWy^v (85) 

71 = -ifl2y™^ 2 (86) 

72 = Shyl^-Ay^ 3 . (87) 

Here h" (0; k,u>, s) tends to as k — > ±00 and uj — > ±00, a property outlined in H3.ll that holds also 
for higher order approximations. In equation (|83[) . the parameter s appears as a factor thanks to the 
homogeneity of the /i-equation (|62p , and affects only the general amplitude of h"(0; k, u), s), which remains 
unknown because of the perturbative approach. Thus, the real and imaginary parts of h"(0;k,u,s)/s 
can be easily plotted after setting a value for y m , and even in this simple model some features of the 
DR map appear to be already present: actually, the imaginary part exhibits a oblique hill and a parallel 
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valley of symmetric heigth separated by a neutral line where h'l(Q) = 0. The equation of this line can be 
deduced from eq. (|83p , and in this simple approximation it turns out to be a linear relation, 



4(i/-y TO 6i) 

k = u> , 

a>iy m v 

where the steepness depends on y m progressively more slowly as y m grows, revealing the dependence on 
the outer boundary location y rn . The hill- valley structure is bounded for very large values of k and uj, 
i.e. this structure vanishes in the regions of the kui plane far from the origin, after eq. At the 

same order in y, i.e. with an approximation for h up to y 5 , but with closure by condition $i{h(y m )} = 0, 
formally expressed as h(y m ) — is where s is a real number, similar results are found, except that the 
roles of real and imaginary part are exchanged, as well as their signs. 

The approximation for h truncated at the next level, i.e. 6th order in y, can be investigated to test 
the dependence of the results on the approximation order. In this case the problem closed by condition 
^{h'(y m )} — gives an expression for h"(0; k, u, s) of the kind 

h"(Q;h,LJ,s)<x — —- — 1 (89) 

y m a + fix oj + p 2 uj 1 + 7i k + -f 2 « + 74 « 

Here the coefficients a, f3\, ...74 are quite cumbersome, but many properties of the previous order remain 
the same, as h"(0; k, uj, s) —> as k — > ±00 and uj — > ±00. As above, s can be factorized without loss of 
generality, and the real and imaginary parts of h"(Q; k, oj, s)/s can be easily plotted. Even in this case, 
the imaginary part h'/(0) exhibits a oblique hill and a parallel valley, but unlike the previous order this 
hill-valley structure is non-symmetric, with an hill height greater than the valley depth. This structure 
vanishes fast outside a k uj domain which is remarkably smaller than in the previous case. Furthermore, 
the transition between the hill and the valley is smoother. It can be shown that the 6th order problem 
with the second closure $t{h(y m )} = still gives similar results. 

A study of the successive approximations at higher orders becomes rapidly complicated, and it is 
probably a slowly converging process, requiring at each order a careful check of the dependence of function 
h on the approximations used for uq and Uq and on the extension of the y-domain. This suggests to 
consider a different approach to the original problem, leading to the numerical method exposed in ij3.2l 
However, the power series study suggests also that some properties have a general validity; for example, 
this is the case of the limit h"(0) — > for large |fc| and \u>\ and of the values /i"(0) = and /i"(0) ^ in 
the origin of the ku> plane. The general shape of h"(0) at the iV-th order of approximation turns out to 
be 

sv N 1 

h"(0; k, u,s)<x -— r (90) 

y m PN[k,LU]y m ) 

where P/v(/c, w; y m ), having y m as a parameter, is a fc, w-polynomial whose order depends on the index 
N. 

The power series approach may serve also as a guide for introducing a reconstruction method of the 
DR map, starting from equation (|52[) rewritten as 

DR r (fc, u) oc i [C r ti'(0;k, w) - dti'(0;k, uj) } 
k 

after omitting the dependence on the parameter s. As said above, in the origin of the kuj plane the real 
and imaginary parts of function h"(0) have different behaviours, namely h'((Q) = and /i"(0) 7^ 0. Thus, 
a simple, proper assumption can be C r = 0, in such a way as to set DR r = at (k,uj) ~ (0,0), when 
the wall doesn't oscillate. Then, the dependence of C on k and to can be expanded in a standard power 
series of the kind 

C = C + Ci k + C 2 uj + C 3 k 2 + C\ k uj + C 5 uj 2 + ■ ■ ■ 

and accounting for the condition (|5ip about the regularity of DR as k — > 0, Co and C 2 must vanish, so 
the series can be rewritten as 

C = k(C 1 +C 3 k + C 4 uj + ■■■). (91) 
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In the practical use, this expansion can be truncated, and the coefficients Cj can be fit on known data, 
satisfying the symmetry (k,u>) <-> (—k, — u), in order to obtain a drag reduction map of the kind 

DR r (fc, u>) oc (C u + C 3i k + C 4t uj + ■■■) h"(0; k, uj) . 
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